sampleSize=1000;
outputSize=100000;
lagMax=3; #Max lag for autocor
Z<-seq(length=outputSize, 0, 0); #output standard gaussian
W<-seq(length=outputSize, 0, 0); #output

#generate observed samples
#X=rexp(sampleSize, rate=1);
#X=runif(sampleSize, min=0, max=1);
X=trace1;
sampleSize=length(X);

#sort X to produce strong autocor
#X=sort(X);

#add two extreme value
X[1]=-100000000000;
X[length(X)]=100000000000;

#map to standard gaussian
F=ecdf(X);
Fx=F(X);
Y<-seq(length=sampleSize-2, 0, 0);
for(i in seq(2, length(X)-1, 1))
{
	Y[i-1]=qnorm(Fx[i], mean=0, sd=1);	
}
#remove first and last element of X
X=X[-length(X)];
X=X[-1];

#map X and Y, fit with a 3rd order polynormial for later Z->W
fit_y_to_x<-lm(X~Y+I(Y^2)+I(Y^3));

#ax^3+bx^2+cx+d
a=fit_y_to_x$coefficient[[4]];
b=fit_y_to_x$coefficient[[3]];
c=fit_y_to_x$coefficient[[2]];
d=fit_y_to_x$coefficient[[1]];
func_y_to_x<-function(x){return (d+c*x+b*x*x+a*x*x*x)};

#code for plot 3rd order polynormial fit
#test_y=seq(-3, 3, length.out=1000)
#test_x=func_y_to_x(test_y)
#plot(Y,X)
#points(test_y, test_x, col="red")

#calculate autocor of the mapped standard gaussian
Ro=acf(Y, lag.max=lagMax, type="correlation", plot=FALSE);

#build covariance matrix (inversely)
Ro_N_N=seq(length=(lagMax+1)*(lagMax+1), 0, 0);
attr(Ro_N_N, "dim")=c(lagMax+1, lagMax+1);
for(i in seq(1,lagMax+1,1))
{
	for(j in seq(1, lagMax+1, 1))
	{
		Ro_N_N[i, j]=Ro$acf[abs(j-i)+1];
	}
}
Ro_1_2=Ro_N_N[c(1:lagMax),c((lagMax+1):(lagMax+1))];
attr(Ro_1_2, "dim")=c(lagMax,1);
Ro_2_1=Ro_N_N[c((lagMax+1):(lagMax+1)), c(1:lagMax)];
attr(Ro_2_1, "dim")=c(1,lagMax);
Ro_2_2=Ro_N_N[c((lagMax+1):(lagMax+1)), c((lagMax+1):(lagMax+1))];
Ro_1_1=Ro_N_N[c(1:lagMax), c(1:lagMax)];

#initialize output
mu_matrix=seq(length=(lagMax),0,0);
for(i in seq(1, lagMax, 1))
{
	Z[i]=Y[i];
	mu_matrix[i]=Y[i];
}
attr(mu_matrix, "dim")=c(lagMax,1);

temp1=Ro_2_1 %*% (solve(Ro_1_1));
new_var=Ro_2_2-(Ro_2_1 %*% (solve(Ro_1_1)) %*% Ro_1_2);
var_error=0;
if(new_var<=0)
{
	var_error=1;
}
new_ro=sqrt(new_var);

#calculate new mu and rio
for(i in seq(lagMax+1, outputSize, 1))
{
	new_mu=temp1 %*% mu_matrix;
	Z[i]=rnorm(1, mean=new_mu, sd=new_ro);
	
	for(j in seq(1, lagMax, 1))
	{
		mu_matrix[j,1]=Z[i-lagMax+j];
	}
}

#inverse back to W (actual output similar to X)
W=func_y_to_x(Z);

#plot cdf
plot(ecdf(X), col="red");
plot(ecdf(W), col="blue", add=TRUE);

#calculate autocor
X_acf=acf(X, lag.max=lagMax, type="correlation", plot=FALSE);
W_acf=acf(W, lag.max=lagMax, type="correlation", plot=FALSE);

Y_acf=acf(Y, lag.max=lagMax, type="correlation", plot=FALSE);
Z_acf=acf(Z, lag.max=lagMax, type="correlation", plot=FALSE);

print(Y_acf);
print(Z_acf);
print(X_acf);
print(W_acf);